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Abstract — Simple and efficient numerical procedures using 
singularity cancellation methods are presented for evaluating 
singular and near-singular potential integrals. Four different 
transformations are compared and the advantages of the Radial- 
angular transform are demonstrated. A method is then 
described for optimizing this integration scheme. 

Index Terms — Numerical integration, Integral equations, 
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I. Introduction 

I T is well known that singular potential integrals involved in 
the integral equations of electromagnetics often require 
special numerical considerations for their evaluation. Until 
recently, the method of singularity subtraction was used 
almost exclusively for the evaluation of singular and near- 
singular terms. However, it becomes unwieldy and prone to 
error for bases of higher order, and it is not conducive to 
object-oriented programming practices, particularly in the 
context of multiple operators. 

To extend the capabilities, accuracy, and maintainability of 
general-purpose codes, the subtraction method is being 
replaced in favor of purely numerical quadrature schemes. 
These schemes employ singularity cancellation methods in 
which a change of variables is chosen such that the Jacobian 
of the transformation cancels the singularity. Recently, a new 
singularity cancellation scheme, the arcsinh transformation, 
was presented for handling l/R singularities [1], For singular 
integrals the method not only has several advantages over 
singularity subtraction methods, but also improves on some 
aspects of other singularity cancellation methods such as Polar 
[2] and Duffy [3] transformations. One drawback of the 
scheme, however, is its inability to efficiently calculate near- 
singular integrals. To this end, we have extended the Duffy 
and polar transformations to handle near-singularities, 
and introduced a new scheme referred to as the Radial-angular 
transformation [4,5]. In this presentation we summarize and 
compare the four transformations for l/R type kernels and 
describe a method for optimizing the Radial-angular 
transformation. 


II. Green’s Functions with l/R singularities 


A. Theory 


In this section the arcsinh, extended Duffy, extended radial 
(or polar) and radial-angular transformations are compared. 
To evaluate potential integrals of the form 

- jkR 

i = f A(r ')— — dD (1) 
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where R = |r -r'| is the distance between an observation point 


at r and source points on triangular domains D , a nearby 
observation point is projected onto a parent triangular 
element, as shown in Fig 1. The parent triangle is then split 
into three subtriangles about the projected observation point. 
The geometry of a typical subtriangle with origin at the 
projected observation point is shown in Fig. 2. 

The general form of the transformed integral over the 
subtriangle is given by 
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where the function h{x,y ) is a product of the basis and 


Green’s potential. Ideally, the transform’s Jacobian J would 
exactly cancel the singular (static) part of the kernel, and the 
subtriangle would map into a rectangular domain such that 
repeated Gauss-Legendre integration of low order may be 
used. A summary of the transformations investigated is given 
in Table 1, followed by a synopsis of their actual properties. 


Observation Point ► # 



projected observation point. 
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TABLE 1: 

SUMMARY OF TRANSFORMATIONS 



TRANSFORMATION 

J (u,v) 

INTEGRATION LIMITS 
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arcsinh : 

• Exact singularity cancellation 

• Rectangular integration domain results only for z = 0 and 
very sensitive to variations in z 

• A single sample point integrates static kernel exactly for 
constant source densities for z = 0 

extended Duffy: 

• Extends Duffy to the case z ^ 0 

• Inexact singularity cancellation for any z 

• A rectangular integration domain results only for z = 0 and 
very sensitive variations in z 

radial: 

• Extends polar form to the case z = 0 

• Exact singularity cancellation 

• Non-rectangular integration domain, but is insensitive to 
variations in z 

radial-angular: 

• Exact singularity cancellation 

• Non-rectangular integration domain forz = 0 , but cosh u 
angular variation in Jacobian is compensated for by v v 

radial limit; insensitive to z = 0 variation 

• Single sample point integrates the static kernel exactly for 
constant source densities for z = 0 
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where {i? I ,R t/ }=||z|, ^z 2 +(/!/sin^) 2 | . For A(r')=l and k = 0, 
Eq. 3 becomes 
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If h 0 d then the kernel in Eq.4 is proportional to hj sin 
In order to smooth the angular dependence we let 
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J ccosh u 

are the weights, and J is the Jacobian of the triangle. Hence, 
the integral can be calculated using repeated Gauss-Legendre 
quadrate with weights w { . . The details of how to evaluate 

potential integrals in simplex coordinates is shown in [1] and, 
therefore, will not be repeated here. 


B. Results 

In this section we compare convergence results for a 
subtriangle with geometry defined as in Fig 3. This geometry 
allows one to limit the degrees of freedom and create 
optimized sampling schemes, which are given in the following 
section. 

For any given subtriangle one determines its maximum 
edge length. As Fig. 2 illustrates, this defines vertices 1 and 2 
of the subtriangle with vertex 3 located between the circular 
arcs in Fig. 3. These arcs have radii equal to the maximum 
edge length and are centered on its end points; they define the 
locus 


For z = 0, the arcsinh and radial-angular approaches are 
nearly identical. The radial-angular transformation presents 
the best compromise for both z = 0 and z ^ 0 and, therefore, 
we will look at this method in more detail. 

From Eq. 1 and Fig. 2 we have 
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Fig. 3. Geometry of a subtriangle showing the Inverse Aspect 
Ratio 


length and are centered on its end points; they define the locus 
of the third vertex whose attached edges are no longer than the 
maximum edge. In order to quantify the subtriangle geometry 
we introduce a quantity called the Inverse Aspect Ratio (IAR), 
which is defined as 


IAR=- 


h 2 h 


i 


(8) 


Where l m:a is the maximum edge length. The IAR is 

normalized such that 0 < IAR < 1 . The maximum IAR occurs 
when vertex 3 is at the intersection of the two circles 
(resulting in an equilateral triangle), and the minimum occurs 
as vertex 3 approaches the maximum edge length. 

The convergence of the potential due to a linear source 
density for a self term and near-self term are shown in Figs. 4 
and 5, respectively. The reference values, shown in Table 2, 
were calculated using the arcsinh and radial-angular methods. 
The subtriangle was subdivided until the two transformations 
agreed to at least 11 significant digits. The maximum edge 
length was chosen to be 1/6 X, and a linear source density was 
defined to be unity at node 1 and zero along the opposite edge. 
Note that for self terms the arcsinh and radial-angular 
transformations have nearly identical convergence behavior. 

For the extended radial and radial-angular methods the 
basis functions contain the terms x=pcost/> and y=ps'm</>, 

where p=V/? 2 -z 2 = ■> jR+z ■> jR-z . For z ^ 0, the convergence 

curves flatten because of the \Ir - z term. In order to 
integrate this term accurately a new quadrature scheme was 
developed based on the family of functions 

F={l,y}R-z,R,RjR-z,..]. (9) 

A partial table for this polynomial-root quadrature scheme is 
shown in Table 3. From Figs. 4 and 5 the radial-angular 
transformation has the best convergence of the four 
transformations studied. However, for z = 0 Gauss-Legendre 
points should be used for the radial sampling, while the 
polynomial-root quadrature scheme should be used for z ^ 0. 


C. Optimization 

Based on the results in the previous section, we can now 
outline a scheme for optimizing the radial-angular 
transformation, the results of which are shown in Table 4. 
Numerical experiments have verified that the largest 
integration error occurs when the third vertex (i.e. singularity) 
lies along one of the bounding circular arcs, the 1/R variation 


in the kernel has the largest variation as measured from the 
vertex to the opposite edge. Hence, if convergence is 
achieved for singularities lying along a circular arc, then 
convergence can be guaranteed for vertices along the same 
line parallel to the largest side, i.e. having the same IAR. 

We define the scheme as optimized which has the fewest 
number of sample points used to achieve convergence given: 
1. The IAR, 2. The distance z from the source plane, 3. Self 
term or near-self term within the plane of the triangle, and 4. 
the desired accuracy chosen by the user. In this work we 
studied the range T4R=[lxlO _8 ,l.o] . We chose a maximum 

edge length of l/6k for the parent triangle. Based on 
experience, we define a near-self term to be an observation 
point that lies within the region Ratio < 1, where 

7?flt/o=|r,|/f max and |r.| is the distance from the observation 

point to the centroid of the parent triangle. Therefore, the 
maximum edge length of a subtriangle can exceed 1/6 X. In 
order to place an upper bound on i max for near-self terms we 
observe that if an observation point is lying on a bounding 



Number of Sample Points 


Fig. 4. Self term convergence for IAR = 0.1. 



Fig. 5. Near-self term convergence for IAR = 0.1 
and d = lxlO' 2 X. 


TABLE 2: 

REFERNCE POTENTIAL VALUES 


D 

REAL 

IMAGINARY 

0.0 

0.015220872777 

-0.0023610287866 

1X1 0' 2 

0.015219642767 

-0.0023610286779 
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TABLE 3: 

RADICAL QUADRATURE POINTS 

N Nodes x, Weights w, 

„ 0.12606123086601956 0.3639172365120473 

_ 0.7139387691339825 0.6360827634879527 

0.0450885041 79695364 0.1 3965395980291434 

3 0.34872938419346483 0.45848221271917206 

0.8306719075452189 0.4018638274779136 

0.019532819681463730 0.06236194190019799 

. 0.17339692801497078 0.25969509521658130 

4 0.522956026924229700 0.4069291 3602039693 

0.88905249698491430 0.27101382686282377 

0.0097091631326854 0.0314958290411988 

0.0927420087986999 0.1478177401386180 

5 0.31 587231 39054450 0.2927739741 663490 

0.6431 824779012770 0.3343492761 944000 

0.921 9651 1 061 291 80 0.1 935631 804594320 


TABLE 4: 

4 SIGNIFICANT DIGIT SAMPLING SCHEMES 
FOR SUBTRIANGLE WITH MAXIMUM EDGE = 5A718 



circular arc and IAR—> 0, the centroid of a triangle will cross 
the edge a distance (■ from a vertex. Hence, the 
maximum edge length achieved for a near-self term is given 
by t =1 ., (Ratio + 2/3) . Choosing Ratio= 1 and 

l =l/6A we obtain I =5/l8/l . 

An additional variable studied for optimization was the 
basis function. For a given observation point, an optimized 
sampling scheme was chosen such that convergence was 
achieved for four different basis functions: a constant function 
and three linear functions. The three linear functions 
represent the bases defined at each of the three vertices. In 
order to justify analyzing each integral independently, as 
opposed to a linear combination, consider the error in 


computing each one as defined by 1 1 - 1 , =6', / . where £ j is 

the relative error in computing the integral. If we define the 
relative error of computing a linear combination of integrals as 

N 


then we have by the triangle inequality 

|^y| - max m . (11) 

Defined in this way, we have the simple result that the relative 
error of a linear combination of N terms, each with error 
| £ f | < £ = max |f ( . | , suffers no loss of significant digits. 

Table 4 shows the optimized sampling schemes of near-self 
terms for 4 significant digits. It was determined that the 
sampling schemes for self and near-self terms were similar 
enough that the near-self schemes could be used for both cases 
without a significant loss in efficiency. Note that the 
optimized sampling schemes are applied to a cylindrical, pill 
box region around the subtriangle. The curved brackets in the 
tables represents an open boundary, while the square bracket 
represents a closed boundary. 

III. Conclusion 

Accurate methods for the purely numerical evaluation of 
singular and near-singular potentials with 1 / R singularities 
have been presented. Of the four transformations studied, the 
radial-angular transformation had the best overall 
convergence. An optimization scheme was developed for this 
transformation based on the observation point location and the 
desired number of significant digits. 
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